Dynamic heterogeneity in a glass forming fluid: susceptibility, structure factor and 

correlation length 
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We investigate the growth of dynamic heterogeneity in a glassy hard-sphere mixture for volume 
fractions up to and including the mode-coupling transition. We use an 80 000 particle system to test a 
new procedure to evaluate a dynamic correlation length £,{t): we determine the ensemble independent 
dynamic susceptibility Xii'^) 8.nd use it to facilitate evaluation of ^(t) from the small wave vector 
behavior of the four-point structure factor. We analyze relations between the a relaxation time Tq,, 
X4('''a)i s-nd ?('''«)• We find that mode-coupling like power laws provide a reasonable description of 
the data over a restricted range of volume fractions, but the power laws' exponents differ from those 
predicted by the inhomogeneous mode-coupling theory. We find ^{tc) ~ In(ra) over the full range 
of volume fractions studied, which is consistent with Adams-Gibbs-type relation. 

PACS numbers: 61.20.Lc,61.20.Ja,64.70.P- 



The search for a growing length scale associated with 
the dramatic slowing down of a supercooled liquid's dy- 
namics is an active area of research. In the last decade, 
a growing dynamic correlation length characterizing dy- 
namic heterogeneity was extensively studied in simula- 
tions [THE] and experiments [5HTT]. and was also investi- 
gated theoretically [T^HTT] . Yet there are still important 
and unresolved issues. 

One popular way to quantify dynamic heterogeneity 
is to focus on the fluctuations of particles' dynamics. 
Since the dynamics is determined by two-point func- 
tions, four-point quantities are introduced to character- 
ize fluctuations of dynamics, the so-called dynamic sus- 
ceptibility Xii^) and associated structure factor S4{q;t). 
Roughly speaking, X4(i) measures the total fluctuations 
of the two-point function characterizing particles' dy- 
namics whereas £'4(5; i) is the Fourier transform of the 
spatially resolved fluctuations. Since the total fluc- 
tuations can formally be obtained by integrating the 
spatially resolved ones, one would naively expect that 
X4(i) = limq^Q Si{q;t). Alas, in a typical simulation at 
least some total fluctuations are not allowed {e.g. the to- 
tal number of particles is fixed); thus Xiit) nieasured in 
simulations is not equal to limg^o 5*4 {q; t) . To distinguish 
the susceptibility determined in an ensemble with quan- 
tity X fixed we henceforth use the symbol X4(i)U- The 
difference between X4(0U and limq_j.o Si{q]t) makes de- 
termination of the dynamic correlation length ^(i) from 
the small wave vector behavior of 5*4(9; ^) more demand- 
ing. It necessitates [5J [3] using significantly larger sys- 
tems than was customary in early simulations. 

In a very interesting development, the difference be- 
tween XiiP) and Xi{t)\x was used by Bcrthicr et al. 
[S] to determine an experimentally accessible bound for 
Xi{t)- Specifically, by using the formalism developed in 
Ref. [18], they showed that XA{t) = X4(i)L + A'(i) where 
X{t) is a correction term due to fluctuations suppressed 
in the constant x ensemble. Importantly, while X4(^)l:r 



cannot be directly determined in experiments, the cor- 
rection term X[t) can. Since X4(*)la; > Oj X{t) provides 
a lower bound to Xi{i)- Both X4(^)l2; and X{t) were cal- 
culated using computer simulations [Ml [19] and it was 
found that X{t) becomes the dominant term close to the 
so-called mode-coupling transition. However, it has not 
been verified that the sum of these two terms agrees with 
the independent extrapolation of S4{q;t) to q = and, 
somewhat surprisingly, the sum has not been used to fa- 
cilitate the evaluation of the dynamic correlation length. 

We note here two difficulties with our present under- 
standing of dynamic heterogeneity. First, there seems to 
be no consensus regarding the scaling relation between 
the length measured at the a relaxation time, f (r^), and 
the relaxation time Ta, even for the range of times ac- 
cessible in computer simulations. Upon approaching the 
mode-coupling transition almost all simulations ^20J find 

1 /z 

a power law ^(r^) ^ Ta ■ However, the range of the scal- 
ing exponents reported is surprisingly large; 1/z varies 
from 0.43 llj to 0.13 [H]. More importantly, it is diffi- 
cult to reconcile relationships between ^(tq) and Ta ex- 
hibited by the simulation results with the experimentally 
determined dynamic correlation lengths p2] . Specifi- 
cally, naive extrapolations of simulational trends result 
in lengths that are orders of magnitude too large for re- 
laxation times at the glass transition temperature |23j . It 
has been shown [TUl [H] that the growth of the dynamic 
susceptibility with the relaxation time slows down near 
the mode-coupling transition. This suggests (but does 
not prove) a similar behavior of the correlation length. 

In this Letter we address the issues mentioned in the 
preceding paragraphs. First, we calculate the ensemble 
independent four-point susceptibility and show explic- 
itly that it agrees very well with the extrapolation of 
54(5; i) to zero wave vector. Next, we analyze the small 
wave- vector behavior of 6*4 ((j;i) and determine £,{t). We 
demonstrate the practical advantage of using X4(^) for 
limq_>.o S4{q; t). Finally, we analyze relations between the 



a relaxation time t^, X4(''a)j ^tnd ^(tq,). Importantly, we 
find the slower-than-power law growth ^(tq,) ~ ln(rQ,). 

We simulated a 50:50 binary mixture of hard spheres 
with Monte Carlo dynamics introduced in Ref. [19]: the 
larger sphere's diameter d-i is 1.4 times larger than the 
smaller sphere diameter d\ , and the dynamics consists of 
random trial moves in a cube of length 0.1 d\. We stud- 
ied systems at fixed numbers of small and large particles 
(A^i and Ni, respectively) or, equivalently, at constant 
volume fraction (\) = (iVidf + N^'^T^jiS^V) iy is the sys- 
tem's volume) and concentration c = Ni/N. We ran four 
trajectories with N = 80 000 particles at = 0.4, 0.45, 
0.5, 0.52, 0.55, 0.56, 0.57, 0.58, and 0.59, and four trajec- 
tories with TV = 10 000 particles at = 0.54, 0.575, 0.58, 
0.585, and 0.59. We ran additional simulations to calcu- 
late derivatives with respect to and c. For (p < 0.585 
our runs are least IOOtq long (ra is defined later) and 
for = 0.59 we ran for 50ra for the 80 000 particle sim- 
ulations and 85 Tq for the 10 000 particle simulations. 
At each (j) at least lOr^ were discarded for equilibration. 
Results are presented in reduced units: time in Monte 
Carlo steps (a Monte Carlo step is one attempted move 
per particle) and lengths in the smaller particle's diame- 
ter. This system was shown (19) to reproduce very well 
the dynamics of an experimental glassy colloidal system. 

To characterize the particles' dynamics we use the 
overlap function w„(t) = 6 {a — |r„(t) — r„(0)|) where 6 
is the Heaviside step function, r„(i) denotes the posi- 
tion of particle n at a time t and a = 0.3. The average 
overlap function Fo{t) = N^^ {J2n ^n(^)) encodes similar 
dynamic information as the self intermediate scattering 
function F,{k;t) = TV"! (Enexp{-«k • [r„(t) - r„(0)]}). 
Here and in the following the brackets (...) denote the 
average over the simulational ensemble in which the vol- 
ume fraction and the composition of the system are fixed. 
The sum in Fo{t) is taken over both the large and small 
particles. We define the a relaxation time Tq, as the time 
at which Fo is equal to 1/e, Fo{Ta) = 1/e. We also evalu- 
ate the mean-square displacement of all the particles and 
the self-diffusion coefficient D. Our results for Tq, and D 
are consistent with those of Brambilla et al. [19] . We ffnd 
that mode-coupling theory-like power laws provide good 
fits to the data for 0.55 < < 0.58 and the power laws' 
exponents are close to those predicted by the theory. Fur- 
thermore, we observe deviations from the power law fits 
for (j) > 0.58. Finally, similarly to Brambilla et al., we 
find that our results for 0.50 < (j) < 0.59 can be well fitted 
by a function showing a stronger divergence at a higher 
volume fraction, bexp[A/{(l)o — 0)^] with <j>o — 0.635. 

It should be emphasized that our simulations extend 
up to and include the mode-coupling transition 0c — 0.59 
whereas previous large scale simulations [5] |6] covered a 
temperature range where mode-coupling power laws pro- 
vide good fits, i.e. a temperature range starting approx- 
imately 15% above the mode-coupling temperature Tc. 

To characterize the heterogeneity of the system's dy- 



namics we examine the dynamic susceptibility Xi{t)\4>,c 
and the four-point structure factor S4{q;t), 



x4{t)u,c = N-U{[w{t)r)-{w{t)y 



(1) 



S,iq- 1) - 7V-1 [{W{q, t)W{-q; t)) - {W{q; t)f) . (2) 

In Eq. (IT]) the subscript 0, c indicates that the suscep- 
tibility is calculated in the simulational ensemble with 
fixed volume fraction and concentration, and W{t) de- 
notes the total overlap at time i, W{t) = J2n^n{t)- In 
Eq. (pi) W{q; t) is the Fourier transform of the spatially 
resolved overlap, W[q]t) = ^^ w„(i) exp[— iq • r„(0)]. 
We expect X4(*)U,c 7^ liniq_i.o 5'4((7;i) and we define the 
ensemble independent dynamic susceptibility, X4(t), as 



Xi{t) = \miSi{q;t). 



(3) 



The difference between Xi{t) and Xi{'^)\<t>,c originates 
from the fluctuations of the volume fraction and the con- 
centration. As recognized by Berthier et al. [^ the con- 
tributions of these fluctuations to Xi (t) can be evaluated 
following Ref. [T5]. We obtain 

X4W « X4WU,c + (yX0W)'gi + ^x0WxcWG2 

+i^o'WG3, (4) 

where p is the number density and Xxit) = dFo(t)/dx. 
There are other terms that contribute to Eq. (|4|, but 
they can be neglected for our system at every studied vol- 
ume fraction. In Eq. (H , G„ are functions of the partial 
structure factors Sapiq) = (NaNp)^-'-/^ {pa{q)pp{—c[)) 
where pa (q) is the Fourier transform of the microscopic 
density of the component a, and a, f3 = 1,2. Explic- 
itly, Gi = dlxiS[l^ + 2dldly^^E^s[l^ + d1x2Sil^ where 



N^/N, and S^°^ 



lim, 



g-J-O 



Sap{q)- 



Shown in Fig. [Tk is the volume fraction dependence 
of the first two terms on the right hand side of Eq. (B 
and the sum of all the terms calculated at t^. We find 
that the xi term becomes the dominant contribution to 
X4{Ta) as increases fM\ . 

To verify that X4(*) calculated from Eq. (|4| agrees with 
limq_j.o 5*4 ((7; t) and to determine ^(t) we need to analyze 
the four-point structure factor. To this end we used the 
80,000 particle simulations to fit Eq. (pi) to several func- 
tions that are based on the following form, 



S4,iq;t) 



A 



C 



i + (e(t)g)2 + i?v [i + (^(t)qy 



(5) 



Specifically, we used (1) the Ornstein-Zernicke (OZ) func- 
tion, i.e. we set i? = and C = 0; (2) a function sug- 
gested by the form of a three-point susceptibility of Ref. 
[13], i.e. we set C = 0; (3) a function suggested by field- 
theoretical considerations of Refs. [T3J[T3JI1S], i.e. we set 
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FIG. 1: (a) The first two terms contributing to X4{'''a) given 
in Eq. Q. Also shown are Xii^a) obtained from Ornstein- 
Zernicke(OZ) extrapolation of §4(5; Tq) to q = and the sum 
of all the terms in Eq. (Ul. (b) The ratio X'i{'''a)/ X'i{'''a)\ ^■ 



A = X4{t)\d, c ^^"^ B = 0; (4) a function utilized in Ref. 
[H], HS4{q;t)] = ln[^] - [^{t)q]'^ + Dq*. All procedures 
resulted in the same lini^_j.o 5*4(9; t) to within error. The 
results for the OZ fits are shown in Fig.[l]as filled squares. 
They agree very well with the open squares which show 
the right-hand-side of Eq. Q . 

Having verified the consistency of Eqs. ^ and Q, we 
now discuss the length ^(tq). Fitting procedures (1) and 
(2) resulted in the same length. Procedure (4) agreed 
with (1) and (2) if we restricted it to wave-vectors q < 
l/S,{Ta)- As expected, procedure (3) resulted in a smaller 
length. This length is approximately 1.2 times smaller 
than the lengths obtained using other fitting procedures, 
independently of (p. Since we established the consistency 
of Eqs. (Is]) and Q, we redid the fits determining £,{Ta) by 
using the right-hand-side of Eq. Q for limg_j.o S4^{q; r^). 
This improved the quality of the fits and reduced the 
uncertainty in ^(tq). As a final revision, since the OZ 
function only provided a good fit for ^{Ta)q < 1.5, we 
restricted the fits to values where q < 1.5/£,{Ta). 

In Fig. [2k we show the results of fitting procedures (1), 
(2), and (4). We also show 1.2^(tq.) when ^ is obtained 
from procedure (3). All these fits produce indistinguish- 
able results. Finally, in Fig. [2J3 we show a scaling plot 
of 5*4 (9; t) shown with the OZ function (solid line). We 
find excellent overlap for every volume fraction and we 
observe deviations from the OZ form only for S^q > 1.5. 
We conclude that using the right-hand-side of Eq. (|4| for 
limq_j.o S4{q;Ta) allows us to reliably determine ^(tq). 

Having tested the procedure to calculate the dynamic 
correlation length we now demonstrate that it can be 
used to get ^(tq) using a moderately large system. We 
find that the correlation lengths determined using 80,000 
and 10,000 particles are virtually identical [see Fig.ls] . 
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FIG. 2: (a) ^(r^) versus (p determined from the different fit- 
ting procedures discussed in the text; the labels denote the 
fitting procedure. The points are statistically the same, and 
the large error bars for (,{Ta) are due to the fits to ln[S'4(g; Tq)]. 
(b) Scaling plot of £'4(9; To)/ Siiq = 0; To,) versus q^{Ta). The 
solid line is the Ornstein-Zernicke function used for the fits. 



Recently, Karmakar et al. O [8] advocated using fi- 
nite size scaling to find the dynamic correlation length 
and showed that C(tq,)'s obtained from this procedure 
are consistent with those determined from the analysis 
of 54 (g; Ta) obtained from very large scale simulations 
(up to 351232 particles). Finite size scaling is attrac- 
tive since it does not require large simulations. However, 
the version used by Karmakar et al. utilizes ensemble- 
dependent quantities, X4:i'''a)\T,n,c {n is the total num- 
ber density) and the fourth central moment of the total 
overlap W(t). Karmakar et al. found that in the tem- 
perature range investigated X4(t4)/ X4(t4)|j,^ ^ « 1.4. 
We believe the temperature independence of this ratio 
might be necessary for the finite size scaling procedure to 
work. Figure [1]d shows that for a range of i^'s the ratio 
X4(''4)/ X4(''4)L c is approximately constant. However, 
upon approaching (j)c it increases rapidly. Thus, Kar- 
makar et a/.'s finite size scaling procedure might not work 
at temperatures close to and below the mode-coupling 
transition Tc. 

Next, we examine scaling relations between t^, X4(''"a) 
and ^('''a)- As shown in Fig. [3k, for ^(r^) > 2 we find 
X4('''q) ^ C('''q)^~'' with 2 — 77 ~ 2.9. Furthermore, as 
shown in Fig. [Sb, for the same range of volume fractions 
where we find good power law fits to r^ and D, we observe 
an approximate power law ^(tq) ~ Ta with 1/z « 0.21. 
The scaling exponents disagree with the inhomogeneous 
mode-coupling predictions |T3j and with Ref. [5]. Inter- 
estingly, 1/z is consistent with some of the earlier studies 
[l m i]. Finally, we find that ^(r^,) ~ In(rc) over the 
whole range of volume fractions studied, thus there is a 
slower-than-power law increase of the dynamic correla- 
tion length with the relaxation time. 

We note that Adam-Gibbs [5S] and Random-First- 
Order- Transition [2 7) theories postulate an exponential 
dependence of the relaxation time, r, on a length, ^, 




FIG. 3: (a) The dynamics susceptibility Xi{'''a) versus the 
correlation length ^(tq). (b) The correlation length ^(r^) 
versus Tq . The straight lines are power law fits over the range 
corresponding to 0.55 < 4> < 0.58. The dashed line is a fit 
to ^ ~ In(Ta) over the whole range of (p. The circles and the 
squares are the results from the 80 000 and the 10 000 particle 
simulations, respectively. 

characterizing the size of correlated regions, r ~ exp(^'''), 
and a relation between ^ and the so-called configurational 
entropy Sc, C ~ (l/5'c)^/(''-^) EH]- A recent study [29j 
confirmed these relations (although somewhat indirectly) 
and found ■0 sa 1, which is consistent with our relation 
between Tq and the dynamic correlation length. This sug- 
gest an intriguing connection between the length charac- 
terizing the size of correlated regions and the dynamic 
heterogeneity length [30] but we leave it for a future 
study. We also acknowledge that another dynamic length 
characterizing single-particle motion was also found to be 
a linear function of \n(Ta) |31| . 

Finally, we examine the (j) dependence of the corre- 
lation length. We find that, in the same range where 
Ta and D are well described by power laws, ^(tq) is 
also well described by a power law with an exponent of 
7^ = 0.46 ±0.03, different from the mode-coupling expo- 
nent of 0.25. The mode-coupling-like fit breaks down 
at the higher volume fractions. Since we found that 
Ta ^ e'^^ and r^ ~ exp{A/((/)o — f/')^}, we fit C('^a) — 
^0 + -^/{4>o ~ i'Y' ■ This function provides a good fit for 
the whole range of cf). It results in ^o = 0.635 ± 0.003. 

To summarize, through large scale computer simula- 
tions we verified that the procedure proposed by Berthier 
et al. results in Xii^) that agrees very well with the in- 
dependent extrapolation limg_j.o Si{q; t). This allowed us 
to propose a new, computationally easier procedure to 
evaluate the dynamic correlation length. Importantly, 
we find a slower-than-power law growth of ^(tq) with t^. 

We thank L. Berthier for discussions and L. Berthier, 
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We gratefully acknowledge the support of NSF Grant No. 
CHE 0909676. 
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